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We study the time evolution of a chain of nonlinear oscillators. We focus on the fractal features 
of the spectral entropy and analyze its characteristic intermediate timescales as a function of the 
nonlinear coupling. A Brownian motion is recognized, with an analytic power-law dependence of its 
diffusion coefficient on the coupling. 
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INTRODUCTION 



II. 



THE SYSTEM 



A system composed of a large number of particles 
(generically) relaxes toward an equilibrium state that is 
independent of the details of the initial state. This is 
one of the fundamental hypotheses of statistical mechan- 
ics. However, from the point of view of Hamiltonian dy- 
namics, the detailed features of the relaxation are not 
thoroughly understood. One of the unsolved fundamen- 
tal questions is how equilibrium is approached when the 
underlying microscopic dynamics is sufficiently chaotic, 
without introducing any randomization or coarsegrain- 
ing "by hand." 

After the pioneering work on the time evolution of non- 
integrable systems pj and ergodicity Q , it is now clear 
that several additional factors play a jprimaryrole in char- 
acterizing the dynamical evolution [3,E1ISl3- However, 
for a large number of particles the KAM argument be- 
comes less effective 0, and the Nekhoroshev bound Q 
for the equilibrium time appears too weak in compari- 
son with numerical results. This scenario has motivated 
a number of numerical studies of Fermi-Pasta-Ulam-likc 
models, in the attempt to clarify the dependence of the 
equilibrium time (defined in terms of suitable indicators) 
on the strength of the nonlinear terms and the number 
of particles y| ■ These studies yield stretched-exponential 
relaxation laws, enforcing the picture that macroscopic 
equilibrium could be built out of local ones Q . 

The aim of this article is to investigate this issue by an- 
alyzing, both analytically and numerically, the dynamics 
of a chain of N coupled anharmonic oscillators at in- 
termediate timescales (for states that are close to equi- 
librium). One of our main results is that at these in- 
termediate timescales the system performs a Brownian 
motion with a diffusion constant that can be accurately 
estimated and turns out to be analytically diverging in 
the coupling constant: as a consequence, a perturbative 
approach to this problem appears sensible. 
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We will study a Hamiltonian made up of an inte- 
grable part and a (small) nonlinear perturbation. This 
is a classical nonintegrable system, that does not possess 
enough integrals of motion. The conjugate variables are 
(q,p) = (<?i, ...qjsr,Pi, ...,Pn), where periodic boundary 
conditions qN+i = Qi, Pn+i = Pi are understood. We 
take N = 2 7 = 128, for convenience of the numerical al- 
gorithm (we observed no significant difference for larger 
N). The Hamiltonian ((/> 4 model) reads 



H(q,p)=H (q,p)+gV(q), (g > 0) 



(1) 
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The quadratic part Hq is easily diagonalized by means of 
a discrete Fourier transform, in terms of the 2A^ normal 
variables q K , p K , with k = (k,a), where k — Q,...,N/2 
and a = 0, 1, a = (a = 1) corresponding to the cosine 
(sine) transform with wave number k. With this coordi- 
nate change Hq becomes 



Hq — E K , E K 
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Pi 
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with the frequency spectrum 
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iin <Uk< WMAX = V 4 + TO 2 . 
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In this article we always set m 2 > 0.1. The value of m 
determines the width of the spectrum and has profound 
consequences on the dynamics, in particular at small non- 
linearities: metastable states, like solitons and breathers, 
are born more easily at small m and this can have drastic 
consequences, both at intermediate and large timescales 
of the equilibration process . 

We integrated the Hamilton equations deriving from 
via a fourth-order Runge-Kutta algorithm in double 



2 



precision. Energy conservation is verified at least up to 1 
part in 10 7 during the whole running time. Such a preci- 
sion is necessary in order to assure that the fluctuations 
due to numerical integrations be negligible with respect 
to the physical ones in which we are interested. 



III. 



SPECTRAL ENTROPY 



From the numerically integrated solutions we analyze 
the behavior of the spectral entropy 



(0) 



where Eo = Ho is the unperturbed total energy. When 
g = all the E K s, and therefore S, are constant. As soon 
as the nonlinearity is switched on, g > 0, the spectral 
entropy becomes a nontrivial function of time. 

The purpose of this letter is to explore the dynamics of 
the system over intermediate timescales. Previous stud- 
ies @ E3 mainly concentrated on the long-time behavior 
of the equilibration process, starting from states that are 
far to equilibrium and smoothing S over sufficiently long 
time intervals. Relaxation from states that are close to 
equilibrium has been more seldom studied In our 

simulation we always set the initial conditions with the 
actions (unperturbed energies) randomly picked from a 
microcanonical ensemble and the angles completely ran- 
dom. Therefore, the fluctuations of S(t) will be studied 
close to equilibrium. 

S displays wild time fluctuations, over a wide range of 
frequencies. We will show that useful, uni vocal informa- 
tion can be obtained from such an irregular function: we 
will first recognize a Brownian structure and then look at 
the characteristic timescales and study their dependence 
on the nonlinear coupling. 

Let us first discuss some analytical properties. At 
equilibrium, average quantities and statistical proper- 
ties should depend only on integrals of motion. For the 
Hamiltonian one argues that the only global integral 
of motion is the total energy E = H or, equivalently, 
the energy per mode e = E/N. The equations for qi 
possess the following scaling symmetry: if qi is a solu- 
tion of the Hamilton equations with coupling g, then 
q[ — yAqi is solution of the Hamilton equations with 
coupling g' = g/A. With this rescaling, the energy is 
changed to E' = AE. A function X, representing the 
average of some quantity at equilibrium and having di- 
mension [X] = length", can only depend on E and g, 
whence 



X{AE,g/A)=A v ' 2 X{E,g), 



(7) 



Therefore, if v = 0, X depends only on the dimcnsionlcss 
product 



that may be considered as the effective strength of the 
nonlinearity. We will study the dynamical properties of 
our system for 1CP 4 < x < 1. 

Average quantities in the weakly nonlinear regime can 
be calculated by using the microcanonical distribution. 
The use of this distribution is motivated by the empir- 
ical observation that the total unperturbed energy Eq, 
although not constant in time, fluctuates less than 1% 
around its mean value (even smaller fluctuations are ob- 
served for very small nonlinearities). We can conclude 
that the primary role of the perturbation V is only to 
allow "collisions" between normal modes (phonon ex- 
change) , provoking in this way the transition to equi- 
librium |17| . without "storing" any significant energy. 
Therefore we can assume E ~ Eq in the following. 

For the spectral entropy at microcanonical equilibrium 
one finds, after a straightforward but somewhat lengthy 
calculation, 

S = <P(N + 1) - (1 - 7 ) = In N - (1 - 7 ) + O (1 V (9) 
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3 + 3N + (tt 2 - 6)iV 2 



3N 2 {N + 1) 
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0.289 

N 



(10) 



where ip is the Euler digamma function [12j and 7 ~ 
0.5772 the Euler-Mascheroni constant. The asymptotic 
behavior of @ is in agreement with previous results 
Eli obtained at canonical equilibrium. In fact, as 
shown in the Appendix, the microcanonical and canoni- 
cal averages of any function of the spectral entropy coin- 
cide. Equations © and (fTUf) are therefore valid also at 
canonical equilibrium. These analytical results are well 
confirmed by numerical simulations and provide a good 
test of the fact that our numerical sample was represen- 
tative of an equilibrium situation and free from "trend" 
components. 



IV. 



CORRELATION FUNCTION 



We study the fractal dimension and the characteristic 
timescales of the entropy by looking at the correlation 
function C for the generalized Brownian process S: 



C(r) 



lim — 



dt (S(t + T) - S(t)Y 



(11) 



x = ge = gE/N, 



(8) 



In general, one can identify a fractal, or generalized 
Brownian motion, by the dependence C cx t 2H . The 
exponent H is related to the fractal dimension by Df — 
2 — H . When using the function one should take 
care of "detrending" S 14]. However, since the system 
starts at equilibrium, where S fluctuates about its con- 
stant mean value © , no detrending is required. We em- 
phasize, however, that we obtained the same results also 



in nonequilibrium situations (not too far from equilib- 
rium), provided the trend component of S(t) was suitably 
removed. 

For a Browni an p rocess one expects a linear depen- 
dence of C on r [Hill, i.e. D f = 3/2, 



C(t) oc t. 



(12) 



Brownian motions are useful idealizations to describe 
physical processes in a simple and coherent mathematical 
way. However, in order to treat an analytic function (like 
S) as a Brownian process, one must identify (at least) one 
timescale, say t\. This timescale is such that by "ob- 
serving" the function at timescales r < t\ one obtains 
a smooth function, while by observing it with a time- 
resolution t > Ti, one recognizes a Brownian process. 
It is possible to unambiguously identify the timescale n 
since C ~ r 2 for sufficiently small r: 



C(r) 



lim — 



1 



* [S(t)T + -S(t)r ' + 0{r 6 ) 



S(t) 



V 2 



0(r 4 ). 



(13) 



At larger r the quadratic dependence changes into the 
linear one Ijl2[l : the timescale at which this change takes 
place is T\. See inset in Figure^ It is important to stress 
that Ti is nothing but the linear timescale for phonons, 
of the order of an inverse characteristic frequency (JSJ of 
the oscillators 
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Up to this timescale, one is able to observe the micro- 
scopic details of the motion in phase space. At larger 
timescales, the motion appears very irregular and the 
microscopic details are lost. 

For a Brownian process which is also bounded another 
timescale appears, since S bounded implies C bounded 
and the law Car must break down at a certain time. 
Let us call this timescale 72 • It is easy to show that 
at equilibrium, for sufficiently large r we have C(r) ~ 
25S 2 , so we can interpret r 2 as the time at which the 
autocorrelation of the entropy vanishes: 

{S(t+T)S(t))-(S(t+T))(S(t)) ~ 0, forr > r 2 . (15) 

This means that S(t + r) and S(t) can be considered 
uncorrelated random numbers chosen from a sample of 
mean S and variance 6S 2 . In terms of motion in phase 
space one can argue that the system starts at t = at (or 
very close to) equilibrium and at t ~ r 2 it has explored a 
sufficiently large part of the equilibrium region, such that 
the microcanonical averages <j9l- Hl(J|) can be used. The 
theoretical predictions (JTJJ [with N = 128], and 
(|13fl are very well verified by the numerical data shown 
in Figure ^ 

The relaxation of the system, when it starts from a 
state far to equilibrium, takes place on a timescale T re i ax 
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FIG. 1: C versus ujt. Notice the quadratic region extending 
up to times of order n (inset) and the linear region in the 
range Ti < T < T2. The saturation level 28 S 2 ~ 0.0044 is in 
full agreement with [for N = 128]. We set x = 0.09 and 
m 2 = 5, so that Q — 2.62. Observe that T2 ~ 200 Ti is clearly 
an intermediate timescale, orders of magnitude smaller than 
Treiax- Larger ratios ts/ti are observed for smaller values of 
x, but the chain of inequalities 11911 always remains valid. 



that can be defined in terms of the effective number of 
excited modes 

n e = exp(S) : (16) 

clearly, if the system is initially far from equilibrium, 

An e (T rc i ax ) ee n e (T rc i ax ) - n e (0) = O(N). (17) 

On the other hand, due to (O-I^UJ, for a system close to 
equilibrium 



An e (r 2 ) 



SS e b =0(VN). (18) 



In this sense, r 2 is an intermediate timescale, character- 
izing, as we have seen, local fluctuations in phase space. 

Finally, if one considers that only a few oscillators can 
exchange energy in a time Ti, one can summarize the 
above discussion by writing 

An e (n) = 0(1) < An e (r 2 ) 

= 0(v / 7V) « An e (T rc i ax ) - 0(N). (19) 



V. DIFFUSION COEFFICIENT AND 
INTERMEDIATE TIMESCALES 

The presence of the linear region (|12l) for the corre- 
lation function C is observed in the whole range of x 
investigated. This enables us to define a diffusion coef- 
ficient D as the rate at which C increases in its linear 
regime, so that in the region n < r < r 2 we have 



C = Dt + C , 



(20) 



where D and Co have in general a nontrivial dependence 
on x. In the following we will perform a systematic study 
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FIG. 2: D versus x — ge — gE/N for different masses. The 
lines are the power law \Y2 lt and the error bars are included 
in the size of the dots. 
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FIG. 3: a = D/x 2 versus m 2 [in the quadratic region 1211 1. 
The solid line is theoretical prediction for large masses a oc 
m -7 . The dashed line is the fit a oc m -3 ' 3 . 



of D, which is the physically most relevant quantity and 
characterizes the Brownian fluctuations. Notice that Cq 
is related to the initial quadratic region jT^ and therefore 
depends on the microscopic details of the motion. 

The diffusion coefficient D has no length dimension 
[v = in J7J)] and therefore one expects it to be only a 
function of x = ge = gE/N (at fixed N and m). This 
expectation is numerically confirmed with high accuracy. 
In particular, D is a monotonically increasing function of 
x, as can be seen from Figure^ Each point in the figure 
is obtained by averaging at least 5 numerical solutions 
with the same value of x. For x < 0.1, D is accurately 
fit by a power law 

D = ax 13 , P= 1.987 ±0.040. (21) 

This is an indication that the intermediate dynamics at 
these timescales can be tackled by a perturbative ap- 
proach. Notice that, if one endeavors to fit the curves in 
Figure [3 with a stretched-exponential law of the type 
D oc exp(const • x~ s ), one finds the very small value 
6 = 9.4 • 10~ 4 . In our opinion, this is a rather strong 
indication in support of a power-law behavior. On the 
other hand, for large x, D saturates to a constant (m- 
independent) value. Observe that the coefficient (3 in 
(|21[) does not depend on to; the m-dependence of a will 
be analyzed in the following (see Figure 0. 

The intermediate timescale r 2 is strictly related to the 
diffusion coefficient D, via the saturation value 2SS 2 in 
l|l()|l. as discussed in connection with Figure A consis- 
tent and natural definition is r 2 = 25S 2 /D, so that for 
x<0.1 

r 2 oc x' p ~ x~ 2 . (22) 

Note that in order to measure r 2 it is not even necessary 
to perform a numerical integration of the Hamilton equa- 
tions for times of order r 2 . It suffices to integrate them 
up to times larger than n (<C t%), get D and hence T2. 

The analytic dependence (|22|l of t 2 on x is in full agree- 
ment with previous results and suggests the validity of a 



perturbative approach 11]. One expects that for suf- 
ficiently small x a sensible fraction of the phase space 
is covered with KAM tori, allowing only ArnoPd diffu- 
sion and yielding a consistent suppression of diffusion 
and a rapid divergence of the macroscopic equilibrium 
time. Related analytical and numerical work hints at a 
non-analytic divergence (such as a stretched exponential) 
of the macroscopic timescales (such as T r ciax) f° r 
x — > 0. The fact that the intermediate timescale here an- 
alyzed diverges with a power law (1221 indicates the local 
presence of a Brownian motion (in the region of param- 
eters studied), so that diffusion should be suppressed on 
macroscopic regions of phase space. 

The dependence of a = D/x 2 on m is not trivial. How- 
ever, in the region of validity of (|21|l and for large masses 
m 2 3> 1, one expects a power-law dependence a oc m~ 7 . 
Indeed, for large m 2 we get from @ e ~ m 2 q 2 , so that, 
at fixed e, q scales like 1/m. Therefore the strength g 
of the quartic potential (and hence x = ge) scales like 
to 4 . Considering that [D] = t^ 1 and the characteristic 
oscillation time (|14f> scales like 1/m, an additional factor 
to is obtained, yielding 

a oc to" 7 , for m 2 > 15 . (23) 

As shown in Figure |31 this prediction is well confirmed 
by our numerical data, provided to 2 > 15. For 0.1 < 
m 2 < 15 we numerically found a oc m~ 3 3 , for which we 
offer no explanation. Moreover, for m — * one expects 
a qualitatively different situation, since WMAx/wmin = 
v4 + to 2 /to — > 00, such a ratio representing the avail- 
able primary resonances between normal modes Ac- 
tually, we observed the formation of metastable states for 
to 2 < 0.05 (not shown in FigureOJ: as a consequence, the 
correlation function showed marked oscillations of defi- 
nite frequency, hindering a consistent definition of a dif- 
fusion coefficient. 
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VI. CONCLUSIONS 

The analytic (power-law) divergence (121'li of the inter- 
mediate timescales describing the Brownian motion on 
the constant energy surface suggests that a perturbative 
app roach, based on a Liouville-Fokker-Planck equation 
[l6| should apply, at least locally, and that an eventually 
non-analytic divergence of the relaxation time from far to 
equilibrium should be ascribed to a nontrivial structure 
of the phase space at larger scales. The dependence of 
the intermediate timescale on the mass and the presence 
of long-lived metastable states indicates that the problem 
is very involved also at the (supposedly simpler) level of 
the local dynamics in phase space. In this context, one 
might speculate that the relaxation from states that are 
not too far to equilibrium depends on the "mesoscopic" 
features of phase space. (The expression "mesoscopic" 
is not used here in its most familiar meaning, related to 
the interplay of classical and quantum effects, but rather 
in the sense of intermediate between "microscopic" and 
"macroscopic," i.e. pertaining to the total system.) The 
approach we propose enables one to extract sensible in- 
formation from the dynamics at intermediate timescales 
(from close to equilibrium) that are usually less studied 
than the longer ones, describing the relaxation from far to 
equilibrium. This has obvious positive spinoffs from the 
perspective of numerical investigations, as it requires less 
computing time. From a more conceptual viewpoint, we 
have proved the existence of a diffusive process (without 
introducing any randomization "by hand"), that consti- 
tutes a building block for the global equilibration process. 



APPENDIX 

Let us show that, in order to compute the average of 
any function of the spectral entropy S, the microcanoni- 
cal and the canonical ensemble are equivalent. This con- 
clusion is a consequence of the following theorem: 

In an ensemble of N noninteracting integrable classi- 
cal systems, where the energies Ei > are homogeneous 
functions of the action variables Ii, the average of func- 
tions of the kind F{Ei/E), where E = Y] Ej, computed 
according to the Boltzmann distribution, is equal to that 
computed according to the microcanonical distribution. 
Moreover, this average is independent of the temperature 
(3 of the canonical ensemble and the energy £ of the mi- 
crocanonical ensemble. 

We prove this theorem by calculating the canonical 
average of a quantity F (the angle variables are trivially 
integrated over since they do not contribute to the ener- 
gies) 



where 



= A. J d * E J{E l )e~^F 



d N Ie 



-f3E 



(A.I) 



(A.2) 



is the canonical partition function and J — \\dlj / dEi\\ 
the Jacobian. Due to our hypotheses, this is a homoge- 
neous function of the Ei's. We obtain 



1 



Ei 



E 



1 



- d£ e"«£ 



J d N x J(£xi)5 (l - x)j F (an) 



if- 



N-l+Na e -/3£ 



x Jd N xJ{ Xi )8 ^l-^aijj F(xi) 

= ^rf3- N{1+a) T(N(l+a)) 

x jd N xJ{ Xi )5 [l-^xA F( Xi ),(A.3) 

where T is the Gamma function, we defined Xi — 
Ei/J2j an( i used the homogeneity of the Jacobian to 



write 



j{tx i )=e a J{x i ), 



(A.4) 



which defines the quantity a. By repeating the same 
derivation with F = 1, one readily shows that 



T(iV(l + a)) 



(A.5) 



where the Z m is the microcanonical partition function 



Z m = J d N I S{£ - E). 



(A.6) 



Incidentally notice that both sides of (|A.5|I are indepen- 
dent of and £. In conclusion 



£ 



N-l + Na 



J d N x J(x % )5 ^1 - F ( x i) 



— [ d N I S(£ -E)F 

%m J 



This proves our assertion. Equations l|§|l- Hl(J|) are readily 
obtained by explicitly calculating the (simpler) canonical 
average of F = S and F — S 2 . 

As a corollary, one sees that if the EiS are random vari- 
ables distributed with e~ Ei , the variables Xi — Ei/ J^j Ej 
are distributed like S(l — J2j x j)- This provides a fast nu- 
merical recipe to generate microcanonically distributed 
variables. 
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